#!/usr/bin/perl

use strict;
use POSIX;
use DBI;

# A tool to get various details about the C1730 database.

# go through the command line options
my %arguments;
for (my $i=0;$i<=$#ARGV;$i++) {
    if ($ARGV[$i] eq "--epoch-list") {
	$arguments{'epoch-list'}=1;
    } elsif ($ARGV[$i] eq "--make-plots") {
	$arguments{'make-plots'}=1;
    } elsif ($ARGV[$i] eq "--source") {
	# a name of a source to make plots for
	$i++;
	push @{$arguments{'sources'}},$ARGV[$i];
    }
}

my $dbh;

my $enable_debug=0;

if ($arguments{'epoch-list'}==1) {
    # open the blazar list from Massaro et al. 2009
    my @ublazars;
    open(BLZ,"massaro_class/bzcat.dat");
    while(<BLZ>) {
	chomp;
	my $line=$_;
	if ($line=~/^(.{4})\s(.{13})\s(.{10})\s(.{11})\s(.{5})\s(.{4})\s(.{3})\s.*$/) {
	    my $tr=&sexa2deci($3)*15;
	    my $td=&sexa2deci($4);
#	    push @ublazars,{ 'name' => $2, 'ra' => $3, 'dec' => $4, 'class' => $7 };
	    push @ublazars,{ 'name' => $2, 'ra' => $tr, 'dec' => $td, 'class' => $7 };
	}
    }
    close(BLZ);	
    my @blazars = sort hash_sorter @ublazars;

    # open the CGRaBS list from Healy et al. 2008
    print "- reading CGRaBS catalogue\n";
    my @ucgrabs;
    open(CGB,"healey_class/table2.dat");
    while(<CGB>) {
	chomp;
	my $line=$_;
	if ($line=~/^(.{10})\s(.{11})\s(.{11})\s(.{7})\s(.{5})\s(.)\s(.)(.{5})\s(.{4})\s(.{6})\s.*$/) {
	    my $tr=&sexa2deci($2)*15;
	    my $td=&sexa2deci($3);
#	    push @ucgrabs,{ 'name' => $1, 'ra' => $2, 'dec' => $3, 'class' => $10 };
	    push @ucgrabs,{ 'name' => $1, 'ra' => $tr, 'dec' => $td, 'class' => $10 };
	}
    }
    close(CGB);
    my @cgrabs = sort hash_sorter @ucgrabs;

    # open the Veron-Cetty & Veron 2010 lists
    print "- reading Veron-Cetty catalogue\n";
    my @uvcagn;
    my @vcfiles=("qso.dat","bllac.dat","agn.dat");
    for (my $i=0;$i<=$#vcfiles;$i++) {
	open(VCA,"veroncetty_class/".$vcfiles[$i]);
	while(<VCA>) {
	    chomp;
	    my $line=$_;
	    if ($line=~/^\s(.)(.{17})\s\s(.{10})\s(.{9})\s(.)\s(.{6})\s(.{4})\s\s(.{6})\s(.{4})\s(.)(.{5})\s(.{4})\s.*$/) {
		my $tr=&sexa2deci($3)*15;
		my $td=&sexa2deci($4);
#		push @uvcagn,{ 'name' => $2, 'ra' => $3, 'dec' => $4, 'class' => $12, 'cat' => $vcfiles[$i] };
		push @uvcagn,{ 'name' => $2, 'ra' => $tr, 'dec' => $td, 'class' => $12, 'cat' => $vcfiles[$i] };
	    }
	}
	close(VCA);
    }
    my @vcagn = sort hash_sorter @uvcagn;

    # open the ARXA list from Flesch 2010
    print "- reading ARXA catalogue\n";
    my @uarxa;
    open(ARXA,"flesch_class/arxa.dat");
    while(<ARXA>) {
	chomp;
	my $line=$_;
	if ($line=~/^(.{2})(.{2})(.{4})(.{3})(.{2})(.{2})\s\s(.{24})\s(.{3})\s.*$/) {
	    my $tr=&sexa2deci($1." ".$2." ".$3)*15;
	    my $td=&sexa2deci($4." ".$5." ".$6);
	    my $tn=$7;
	    $tn=~s/^\s*//;
	    $tn=~s/\s*$//;
#	    push @uarxa,{ 'name' => $7, 'ra' => $1." ".$2." ".$3, 'dec' => $4." ".$5." ".$6, 'class' => $8 };
	    push @uarxa,{ 'name' => $tn, 'ra' => $tr, 'dec' => $td, 'class' => $8 };
	}
    }
    close(ARXA);
    my @arxa = sort hash_sorter @uarxa;

    # open the 1LAC list from Abdo et al. 2010
    print "- reading 1LAC catalogue\n";
    my @ulac;
    open(LAC,"1lac/table1.dat");
    while(<LAC>) {
	chomp;
	my $line=$_;
	if ($line=~/^.{5}(.{12})\s(.{30})\s(.{9})\s(.{9})\s(.{5})\s(.{5})\s(.{4})\s(.{4})\s(.{3})\s(.{8})\s(.{4})\s(.{3})\s(.{4})\s(.{4})\s(.{5})\s(.)\s(.)\s.*$/) {
	    push @ulac, { 'name' => $1, 'ra' => $3, 'dec' => $4, 'class' => $8, 'sed' => $9, 'z' => $10,
			  'gflux' => $11, 'gfluxerror' => $12, 'gamma' => $13, 'gammaerror' => $14,
			  'variable' => $16, 'speccurve' => $17 };
	}
    }
    close(LAC);
    my @lac = sort hash_sorter @ulac;

    # we query the database to make a list of C1730 sources vs the
    # epochs that they were observed in
    &connectdb();

    my %missing_epochs = ( 
			   '0208-512' => { 'epochs' => [ '2008-10-13 00:00:00' ]},
			   '0414-189' => { 'epochs' => [ '2008-10-13 00:00:00' ]},
			   '0454-234' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '0454-463' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '0537-441' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '0616-116' => { 'epochs' => [ '2008-10-13 00:00:00' ]},
			   '1334-127' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '1424-418' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '1622-297' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '1622-253' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '1730-130' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '1908-201' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '1933-400' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '1934-638' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '1936-155' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '2052-474' => { 'epochs' => [ '2008-10-13 00:00:00' ] },
			   '2155-304' => { 'epochs' => [ '2008-10-13 00:00:00' ] }
			   );			   

    print "- querying database for C1730 information\n";
    my $epoch_query="SELECT fluxdata.source,".
	"to_array(date_trunc('day',time)) FROM fluxdata INNER JOIN ".
	"observation ON (observation.id=ANY(fluxdata.observation_ids)) ".
	"WHERE project_code='C1730' GROUP BY fluxdata.source ORDER BY ".
	"fluxdata.source";
    my $epoch_handle=$dbh->prepare($epoch_query);
    $epoch_handle->execute;
    my %epoch_data;
    while (my @epoch_response=$epoch_handle->fetchrow()) {
	$epoch_data{$epoch_response[0]} = {
	    'epochs' => $epoch_response[1]
	};
	&unique_array($epoch_data{$epoch_response[0]}->{'epochs'});
    }
    $epoch_handle->finish;

    foreach my $k (keys %missing_epochs) {
	if (!$epoch_data{$k}) {
	    $epoch_data{$k} = $missing_epochs{$k};
	} else {
	    push @{$epoch_data{$k}->{'epochs'}},@{$missing_epochs{$k}->{'epochs'}};
	}
    }

    my @all_epochs;
    my $maxnamelength=0;

    # get some RAs and Decs from the database
    foreach my $k (keys %epoch_data) {
	my @rd = &get_radec($k);
	$epoch_data{$k}->{'ra'}=$rd[0];
	$epoch_data{$k}->{'dec'}=$rd[1];
	for (my $i=0;$i<=$#{$epoch_data{$k}->{'epochs'}};$i++){
	    my $e = ${$epoch_data{$k}->{'epochs'}}[$i];
	    $e =~ s/\s+00\:00\:00$//;
	    push @all_epochs, $e;
	    ${$epoch_data{$k}->{'epochs'}}[$i] = $e;
	}
	my $l = length $k;
	if ($l > $maxnamelength) {
	    $maxnamelength = $l;
	}
    }
    &unique_array(\@all_epochs);
    my @sorted_epochs = sort(@all_epochs);
    for (my $i=0;$i<$#sorted_epochs;$i++) {
	my $s = &within_day($sorted_epochs[$i], $sorted_epochs[$i+1]);
	if ($s == 1) {
	    splice @sorted_epochs,$i+1,1;
	}
    }

    open(JSON,">c1730_epochs_radec.json");
    open(TABLE,">c1730_source_latex_table.txt");
    print JSON "{";

    printf "%".$maxnamelength."s ", " ";
    print JSON " 'epochs': [";
    for (my $i=0;$i<=$#sorted_epochs;$i++) {
	print $sorted_epochs[$i]." ";
	if ($i>0) {
	    print JSON ",";
	}
	print JSON "'".$sorted_epochs[$i]."'";
    }
    print "\n";
    print JSON "] ";

    my @sources = keys %epoch_data;
    my @sorted_sources = sort @sources;
    my @table_sources;
    my @json_sources;

    print "- matching fluxes and catalogue information\n";
    for (my $j=0;$j<=$#sorted_sources;$j++) {
	my $k = $sorted_sources[$j];
	# get the fluxes
	my $flux_query="SELECT fluxdata.source,fluxdata.frequency,".
	    "tripleamp_value,tripleamp_scatter_measured,date_trunc('day',time)".
	    " FROM fluxdata INNER JOIN observation ON ".
	    "(observation.id=fluxdata.observation_ids[1]) WHERE ".
	    "fluxdata.source='".$k."' AND project_code='C1730'";
	my $flux_handle=$dbh->prepare($flux_query);
	$flux_handle->execute;
	my %tfluxes;
	while (my @flux_response=$flux_handle->fetchrow()) {
	    # find the nearest epoch
	    my $nep='';
	    for (my $i=0;$i<=$#sorted_epochs;$i++) {
		if (&within_day($flux_response[4],$sorted_epochs[$i]) == 1) {
		    $nep=$sorted_epochs[$i];
		    last;
		}
	    }
	    if ($nep ne '') {
		if (!exists $tfluxes{$nep}) {
		    $tfluxes{$nep}={ 'frequency' => [],
				     'flux' => [],
				     'error' => [],
				     'n' => [] };
		}
		my $fcf=0;
		for (my $cc=0;$cc<=$#{$tfluxes{$nep}->{'frequency'}};$cc++) {
		    if (${$tfluxes{$nep}->{'frequency'}}[$cc]==$flux_response[1]) {
			# the same frequency
			$fcf=1;
			${$tfluxes{$nep}->{'flux'}}[$cc]+=${$flux_response[2]}[0];
			${$tfluxes{$nep}->{'n'}}[$cc]++;
			${$tfluxes{$nep}->{'error'}}[$cc]=(${$tfluxes{$nep}->{'error'}}[$cc]>
							   ${$flux_response[3]}[0]) ?
							   ${$tfluxes{$nep}->{'error'}}[$cc] :
							   ${$flux_response[3]}[0];
		    }
		}
		if ($fcf==0) {
		    push @{$tfluxes{$nep}->{'frequency'}},$flux_response[1];
		    push @{$tfluxes{$nep}->{'flux'}},${$flux_response[2]}[0];
		    push @{$tfluxes{$nep}->{'error'}},${$flux_response[3]}[0];
		    push @{$tfluxes{$nep}->{'n'}},1;
		}
	    }
	}
	$flux_handle->finish;

	push @json_sources, { 'name' => $sorted_sources[$j],
			      'ra' => $epoch_data{$k}->{'ra'},
			      'dec' => $epoch_data{$k}->{'dec'},
			      'epochdata' => [], 'sameAs' => '' };
	my $cjhandle=$json_sources[$#json_sources];
	my ($rarep, $decrep);
	printf "%".$maxnamelength."s ", $k;
	my $uname = uc($k);

	# check if this is the same position as another source
	my $nsource=-1;
	my $h1n=&sexa2deci($epoch_data{$k}->{'ra'})*15;
	my $d1n=&sexa2deci($epoch_data{$k}->{'dec'});
	for (my $i=0;$i<=$#table_sources;$i++) {
	    my $h2n=&sexa2deci($table_sources[$i]->{'ra'})*15;
	    my $d2n=&sexa2deci($table_sources[$i]->{'dec'});
	    my $ad=&ang_dist($h1n,$d1n,$h2n,$d2n);
	    if ($ad<=10/3600) { # separation leq 10 arcseconds
		$nsource=$i;
		if ($table_sources[$i]->{'altname'} ne "") {
		    $table_sources[$i]->{'altname'}.=",";
		}
		$table_sources[$i]->{'altname'}.=$uname;
		$cjhandle->{'sameAs'}=$table_sources[$i]->{'name'};
	    }
	}
	if ($nsource==-1) {
	    push @table_sources,{ 'name' => $uname, 'ra' => $epoch_data{$k}->{'ra'},
				  'dec' => $epoch_data{$k}->{'dec'},
				  'epochs' => [], 'altname' => "", 'epochnames' => [],
				  'class' => "" };
	    # see if this source is classified by Massaro
	    my $mi=&find_val(\@blazars,$h1n,$d1n);
	    if ($mi>-1) {
		$table_sources[$#table_sources]->{'class'}=$blazars[$mi]->{'class'}."\$^1\$";
		$cjhandle->{'bzcat_class'}=$blazars[$mi]->{'class'};
		$cjhandle->{'bzcat_name'}=$blazars[$mi]->{'name'};
	    }
	    # see if this source is classified by Healey
	    $mi=&find_val(\@cgrabs,$h1n,$d1n);
	    if ($mi>-1) {
		$cgrabs[$mi]->{'class'}=~s/\s//g;
		if ($cgrabs[$mi]->{'class'} ne "") {
		    if ($table_sources[$#table_sources]->{'class'} ne "") {
			$table_sources[$#table_sources]->{'class'}.=",";
		    }
		    $table_sources[$#table_sources]->{'class'}.=$cgrabs[$mi]->{'class'}."\$^2\$";
		    $cjhandle->{'cgrabs_class'}=$cgrabs[$mi]->{'class'};
		    $cjhandle->{'cgrabs_name'}=$cgrabs[$mi]->{'name'};
		}
	    }
	    # see if this source is classified by Veron-Cetty
	    $mi=&find_val(\@vcagn,$h1n,$d1n);
	    if ($mi>-1) {
		$vcagn[$mi]->{'class'}=~s/\s//g;
		if ($vcagn[$mi]->{'class'} ne "") {
		    if ($table_sources[$#table_sources]->{'class'} ne "") {
			$table_sources[$#table_sources]->{'class'}.=",";
		    }
		    $table_sources[$#table_sources]->{'class'}.=$vcagn[$mi]->{'class'}."\$^3\$";
		    $cjhandle->{'vcagn_class'}=$vcagn[$mi]->{'class'};
		    $cjhandle->{'vcagn_name'}=$vcagn[$mi]->{'name'};
		}
	    }
	    # see if this source is classified by Flesch
	    $mi=&find_val(\@arxa,$h1n,$d1n);
	    if ($mi>-1) {
		$arxa[$mi]->{'class'}=~s/\s//g;
		if ($arxa[$mi]->{'class'} ne "") {
		    if ($table_sources[$#table_sources]->{'class'} ne "") {
			$table_sources[$#table_sources]->{'class'}.=",";
		    }
		    $table_sources[$#table_sources]->{'class'}.=$arxa[$mi]->{'class'}."\$^4\$";
		    $cjhandle->{'arxa_class'}=$arxa[$mi]->{'class'};
		    $cjhandle->{'arxa_name'}=$arxa[$mi]->{'name'};
		}
	    }
	    # see if this source is in the 1LAC
	    $mi=&find_val(\@lac,$h1n,$d1n);
	    if ($mi>-1) {
		$lac[$mi]->{'class'}=~s/\s//g;
		if ($lac[$mi]->{'class'} ne "") {
		    if ($table_sources[$#table_sources]->{'class'} ne "") {
			$table_sources[$#table_sources]->{'class'}.=",";
		    }
		    $table_sources[$#table_sources]->{'class'}.=$lac[$mi]->{'class'}."\$^4\$";
		    $cjhandle->{'lac_class'}=$lac[$mi]->{'class'};
		    $cjhandle->{'lac_name'}=$lac[$mi]->{'name'};
		    $cjhandle->{'lac_sed'}=$lac[$mi]->{'sed'};
		    $cjhandle->{'lac_z'}=$lac[$mi]->{'z'};
		    $cjhandle->{'lac_gflux'}=$lac[$mi]->{'gflux'};
		    $cjhandle->{'lac_gfluxerror'}=$lac[$mi]->{'gfluxerror'};
		    $cjhandle->{'lac_gamma'}=$lac[$mi]->{'gamma'};
		    $cjhandle->{'lac_gammaerror'}=$lac[$mi]->{'gammaerror'};
		    $cjhandle->{'lac_variable'}=$lac[$mi]->{'variable'};
		    $cjhandle->{'lac_speccurve'}=$lac[$mi]->{'speccurve'};
		}
	    }
	}

	my $tablen=0;
	for (my $i=0;$i<=$#sorted_epochs;$i++) {
	    my $L = length $sorted_epochs[$i];
	    my $n = 1;
	    for (my $l=0;$l<=$#{$epoch_data{$k}->{'epochs'}};$l++) {
		if (&within_day(${$epoch_data{$k}->{'epochs'}}[$l],
				$sorted_epochs[$i]) == 1) {
		    $n = 0;
		    last;
		}
	    }
	    printf "%".$L."s ", ($n==1) ? "no" : "yes";
	    if ($n==0) {
		if ($nsource==-1){
		    push @{$table_sources[$#table_sources]->{'epochs'}}, ($i+1);
		} else {
		    my $nepoch=1;
		    for (my $p=0;$p<=$#{$table_sources[$nsource]->{'epochs'}};$p++) {
			if (${$table_sources[$nsource]->{'epochs'}}[$p]==($i+1)) {
			    $nepoch=0;
			    last;
			}
		    }
		    if ($nepoch==1) {
			push @{$table_sources[$nsource]->{'epochs'}},($i+1);
		    }
		}
	    }
	    my $tjh = { 'epoch' => $sorted_epochs[$i], 'fluxes' => [] };
	    $tjh->{'observed'} = ($n==1) ? 'false' : 'true';
	    if (exists $tfluxes{$sorted_epochs[$i]}) {
		for (my $j=0;$j<=$#{$tfluxes{$sorted_epochs[$i]}->{'frequency'}};$j++) {
		    my $tjf = { 'frequency' => ${$tfluxes{$sorted_epochs[$i]}->{'frequency'}}[$j],
				'flux_density' => ${$tfluxes{$sorted_epochs[$i]}->{'flux'}}[$j]/
				    ${$tfluxes{$sorted_epochs[$i]}->{'n'}}[$j],
				'uncertainty' => ${$tfluxes{$sorted_epochs[$i]}->{'error'}}[$j] };
		    push @{$tjh->{'fluxes'}},$tjf;
		}
	    }
	    push @{$cjhandle->{'epochdata'}},$tjh;
	}
	print "\n";
    }

    print "- writing out LaTeX table for paper\n";
    # sort the table sources by RA
    my @stable_sources=sort table_sorter @table_sources;
    my $first_n=55;
    my $other_n=64;
    my $n_other;
    for (my $i=0;$i<=$#stable_sources;$i++) {
	if ($stable_sources[$i]->{'name'} eq "1934-638" ||
	    $stable_sources[$i]->{'name'} eq "URANUS") {
	    next;
	}
	print TABLE $stable_sources[$i]->{'name'};
#	if ($stable_sources[$i]->{'altname'} ne "") {
#	    print TABLE " (".$stable_sources[$i]->{'altname'}.")";
#	}
	print TABLE " & ".
	    $stable_sources[$i]->{'ra'}." & ".
	    $stable_sources[$i]->{'dec'}." & ";
	my @sepochs = sort {$a <=> $b} @{$stable_sources[$i]->{'epochs'}};
	for (my $k=0;$k<=$#sepochs;$k++) {
	    if ($k>0) {
		print TABLE ",";
	    }
	    print TABLE $sepochs[$k];
	}
	print TABLE " & ".$stable_sources[$i]->{'class'};
	print TABLE "\\\\\n";
	if ($first_n>0) {
	    $first_n--;
	    if ($first_n==0) {
		print TABLE "\n\n";
		$n_other=$other_n;
	    }
	} else {
	    $n_other--;
	    if ($n_other==0) {
		print TABLE "\n\n";
		$n_other=$other_n;
	    }
	}
    }

    print "- writing out JSON file for web and analysis\n";
    # make the JSON output
    for (my $i=0;$i<=$#json_sources;$i++) {
	print JSON ",\n'".$json_sources[$i]->{'name'}."': {".
	    " 'ra': '".$json_sources[$i]->{'ra'}."', ".
	    "'dec': '".$json_sources[$i]->{'dec'}."', ";
	if ($json_sources[$i]->{'sameAs'} ne '') {
	    print JSON "'sameAs': '".$json_sources[$i]->{'sameAs'}."', ";
	}
	my @refcats=("bzcat","cgrabs","vcagn","arxa","lac");
	for (my $j=0;$j<=$#refcats;$j++) {
	    if ($json_sources[$i]->{$refcats[$j]."_class"}) {
		print JSON "'".$refcats[$j]."Class': '".
		    $json_sources[$i]->{$refcats[$j]."_class"}."', ".
		    "'".$refcats[$j]."Name': '".
		    $json_sources[$i]->{$refcats[$j]."_name"}."', ";
	    }
	    if ($json_sources[$i]->{$refcats[$j]."_sed"}) {
		print JSON "'".$refcats[$j]."SED': '".
		    $json_sources[$i]->{$refcats[$j]."_sed"}."', ".
		    "'".$refcats[$j]."Z': ".
		    $json_sources[$i]->{$refcats[$j]."_z"}.", ".
		    "'".$refcats[$j]."GFlux': ".
		    $json_sources[$i]->{$refcats[$j]."_gflux"}.", ".
		    "'".$refcats[$j]."GFluxError': ".
		    $json_sources[$i]->{$refcats[$j]."_gfluxerror"}.", ".
		    "'".$refcats[$j]."Gamma': ".
		    $json_sources[$i]->{$refcats[$j]."_gamma"}.", ".
		    "'".$refcats[$j]."GammaError': ".
		    $json_sources[$i]->{$refcats[$j]."_gammaerror"}.", ".
		    "'".$refcats[$j]."Variable': '".
		    $json_sources[$i]->{$refcats[$j]."_variable"}."', ".
		    "'".$refcats[$j]."SpecCurve': '".
		    $json_sources[$i]->{$refcats[$j]."_speccurve"}."', ";
	    }
	}
	print JSON "\n'data': [";
	for (my $j=0;$j<=$#{$json_sources[$i]->{'epochdata'}};$j++) {
	    my $tref=$json_sources[$i]->{'epochdata'};
	    if ($j>0) {
		print JSON ",";
	    }
	    print JSON "\n{ 'epoch': '".${$tref}[$j]->{'epoch'}."',".
		"'observed': ".${$tref}[$j]->{'observed'};
	    my $fref=${$tref}[$j]->{'fluxes'};
	    if ($#{$fref}>-1) {
		print JSON ", 'observations': [";
		for (my $k=0;$k<=$#{$fref};$k++) {
		    if ($k>0) {
			print JSON ",";
		    }
		    print JSON "\n{ 'frequency': ".${$fref}[$k]->{'frequency'}.
			", 'flux_density': ".${$fref}[$k]->{'flux_density'}.
			", 'uncertainty': ".${$fref}[$k]->{'uncertainty'}." }";
		}
		print JSON "]";
	    }
	    print JSON "}";
	}
	print JSON "]}";
    }
    print JSON "}\n";
}


if ($arguments{'make-plots'}==1) {
    &connectdb();

    # do we have a list of specific sources to make plots for?
    my $plot_query = "SELECT DISTINCT ON (fluxdata.id) fluxdata.id,fluxdata.source,".
	"fluxdata.frequency,tripleamp_value,".
	"tripleamp_scatter_measured,date_trunc('day',time) FROM fluxdata INNER JOIN observation ON ".
	"(observation.id=ANY(observation_ids)) WHERE project_code='C1730'";
    if (!$arguments{'sources'}) {
	# no, so we do all of them

    } else {
	$plot_query .= " AND (";
	for (my $i=0;$i<=$#{$arguments{'sources'}};$i++) {
	    if ($i>0) {
		$plot_query .= " OR ";
	    }
	    $plot_query .= "fluxdata.source='".${$arguments{'sources'}}[$i]."'";
	}
	$plot_query .= ")";
    }
    my $plot_handle=$dbh->prepare($plot_query);
    $plot_handle->execute;

    my %plot_data;
    while (my @plot_response=$plot_handle->fetchrow()) {
	# check for this source
	if (!$plot_data{$plot_response[1]}) {
	    $plot_data{$plot_response[1]} = {};
	}

	# check for this epoch
	my $cday = $plot_response[5];
	$cday =~ s/\s+00\:00\:00$//;
	foreach my $k (keys %{$plot_data{$plot_response[1]}}) {
	    my $s = &within_day($cday, $k);
	    if ($s == 1) {
		# already got one within a day, use it instead
		$cday = $s;
	    }
	}

	# check whether this source has already been registered on this day
	if (!$plot_data{$plot_response[1]}->{$cday}) {
	    $plot_data{$plot_response[1]}->{$cday} = {
		'frequency' => [],
		'flux' => [],
		'fluxunc' => [],
		'n' => []
	    };
	} else {
	    # we have a repeat observation, is it just at a different frequency?
	    my $found_freq=0;
	    for (my $k=0;$k<=$#{$plot_data{$plot_response[1]}->{$cday}->{'frequency'}};$k++) {
		if (${$plot_data{$plot_response[1]}->{$cday}->{'frequency'}}[$k] ==
		    $plot_response[2]) {
		    $found_freq=1;
		    # take the average flux and the biggest uncertainty
		    my $nflux = ${$plot_data{$plot_response[1]}->{$cday}->{'flux'}}[$k] *
			${$plot_data{$plot_response[1]}->{$cday}->{'n'}}[$k];
		    $nflux += ${$plot_response[3]}[0];
		    ${$plot_data{$plot_response[1]}->{$cday}->{'n'}}[$k]++;
		    ${$plot_data{$plot_response[1]}->{$cday}->{'flux'}}[$k] += $nflux/
			${$plot_data{$plot_response[1]}->{$cday}->{'n'}}[$k];
		    ${$plot_data{$plot_response[1]}->{$cday}->{'fluxunc'}}[$k] =
			(${$plot_data{$plot_response[1]}->{$cday}->{'fluxunc'}}[$k] >
			 ${$plot_response[4]}[0]) ?
			 ${$plot_data{$plot_response[1]}->{$cday}->{'fluxunc'}}[$k] :
			 ${$plot_response[4]}[0];
		    last;
		}
	    }
	    if ($found_freq == 0) {
		# add a new flux measurement
		push @{$plot_data{$plot_response[1]}->{$cday}->{'frequency'}},
		$plot_response[2];
		push @{$plot_data{$plot_response[1]}->{$cday}->{'flux'}},
		${$plot_response[3]}[0];
		push @{$plot_data{$plot_response[1]}->{$cday}->{'fluxunc'}},
		${$plot_response[4]}[0];
		push @{$plot_data{$plot_response[1]}->{$cday}->{'n'}}, 1;
	    }
	}
    }

    
}

sub within_day {
    my $date1 = shift;
    my $date2 = shift;
    
    my ($y1,$m1,$d1) = split(/\-/, $date1);
    my ($y2,$m2,$d2) = split(/\-/, $date2);

    if ( $y1 == $y2 && $m1 == $m2 ) {
	if (abs($d2 - $d1) < 2) {
	    return 1;
	} else {
	    return 0;
	}
    } else {
	return 0;
    }
}

sub unique_array {
    my $array_ref = shift;

    for (my $i=0;$i<$#{$array_ref};$i++) {
	my $v = ${$array_ref}[$i];
	if ($enable_debug==1) {
	    print " found value $v\n";
	}
	for (my $j=$i+1;$j<=$#{$array_ref};$j++) {
	    my $u = ${$array_ref}[$j];
	    if ($enable_debug==1) {
		print "  found second value $u: ";
	    }
	    if ($v eq $u) {
		# drop this element
		splice @{$array_ref},$j,1;
		$j--;
		if ($enable_debug==1) {
		    print "dropping array element\n";
		}
	    } elsif ($enable_debug==1) {
		print "not dropping array element\n";
	    }
	}
    }
}

sub connectdb {
    # database details
    my $dsn="dbi:Pg:dbname=atca;host=mentok;port=5432";
    my $user="postgres";
    my $password="";
    
    # Setup the database connection
    if (!$dbh) {
	$dbh=DBI->connect($dsn,$user,$password,{ RaiseError => 0, AutoCommit => 1});
    }
}

sub get_radec {
    my $sname = shift;

    my $radec_query="SELECT ra,dec FROM observation WHERE (source='".$sname."'".
    	" OR source='".uc($sname)."')".
	" AND project_code='C1730' LIMIT 1";
    my $radec_handle=$dbh->prepare($radec_query);
    $radec_handle->execute;
    my @radec_response=$radec_handle->fetchrow();
    $radec_handle->finish;

    my $nra=&deci2sexa($radec_response[0]);
    my $ndec=&deci2sexa($radec_response[1]);
    return ($nra,$ndec);
}

sub sexa2deci {
    my $sexa_number=shift;
    
    my $sign=1;
    if ($sexa_number=~/^\-/) {
	$sign=-1;
	$sexa_number=~s/^\-//;
    }
    my @sexaels=split(/[\s\:]/, $sexa_number);
    my $deci_number=$sexaels[0]+$sexaels[1]/60+$sexaels[2]/3600;
    $deci_number*=$sign;

    return $deci_number;
}

sub deci2sexa {
    my $deci_number=shift;

    my $sign=1;
    if ($deci_number<0) {
	$sign=-1;
    }
    $deci_number*=$sign;

    my $fp=floor($deci_number);
    $deci_number-=$fp;
    $deci_number*=60;
    my $sp=floor($deci_number);
    $deci_number-=$sp;
    $deci_number*=60;
    my $tp=$deci_number;
    
    $fp*=$sign;
    my $sexastring=sprintf "%02d:%02d:%04.1f",
    $fp,$sp,$tp;

    return $sexastring;

}

sub polar2rec {
    my $r=shift;
    my $theta=shift;

    my $x=$r*cos($theta);
    my $y=$r*sin($theta);

    return ($x,$y);
}

sub rec2polar {
    my $x=shift;
    my $y=shift;

    my $r=sqrt($x*$x+$y*$y);
    my $theta=atan2($y,$x);

    return ($r,$theta);
}

sub ang_dist {
    my $ra1=shift;
    my $dec1=shift;
    my $ra2=shift;
    my $dec2=shift;
    my $dtor=3.141592654/180.0;

    my ($rxy,$z1,$x1,$y1,$z2,$x2,$y2);
    ($rxy,$z1)=&polar2rec(1.0,$dec1*$dtor);

    ($x1,$y1)=&polar2rec($rxy,$ra1*$dtor);

    ($rxy,$z2)=&polar2rec(1.0,$dec2*$dtor);

    ($x2,$y2)=&polar2rec($rxy,$ra2*$dtor);

    my $cs=$x1*$x2+$y1*$y2+$z1*$z2;
    my $xc=$y1*$z2-$z1*$y2;
    my $yc=$z1*$x2-$x1*$z2;
    my $zc=$x1*$y2-$y1*$x2;
    my $sn=sqrt($xc*$xc+$yc*$yc+$zc*$zc);
    my ($r,$pDist)=&rec2polar($cs,$sn);

    return ($pDist/$dtor);
}

sub table_sorter {
    my $ar=&sexa2deci($a->{'ra'})*15;
    my $br=&sexa2deci($b->{'ra'})*15;

    if ($ar == $br) {
	my $ad=&sexa2deci($a->{'dec'});
	my $bd=&sexa2deci($b->{'dec'});

	if ($ad < $bd) {
	    return -1;
	} elsif ($ad == $bd) {
	    return 0;
	} elsif ($ad > $bd) {
	    return 1;
	}
    } elsif ($ar < $br) {
	return -1;
    } elsif ($ar > $br) {
	return 1;
    }
}

sub hash_sorter {
    my $ar=$a->{'ra'};
    my $br=$b->{'ra'};

    if ($ar == $br) {
	my $ad=$a->{'dec'};
	my $bd=$b->{'dec'};

	if ($ad < $bd) {
	    return -1;
	} elsif ($ad == $bd) {
	    return 0;
	} elsif ($ad > $bd) {
	    return 1;
	}
    } elsif ($ar < $br) {
	return -1;
    } elsif ($ar > $br) {
	return 1;
    }
}

sub find_val {
    my $listref=shift;
    my $valfind_ra=shift;
    my $valfind_dec=shift;

    my $dvalfind_ra=&sexa2deci($valfind_ra);
    my $dvalfind_dec=&sexa2deci($valfind_dec);
    my ($lowlim,$highlim)=&bin_sort($listref,$dvalfind_ra);

    for (my $i=$lowlim;$i<=$highlim;$i++) {
	my $ad=&ang_dist($dvalfind_ra,$dvalfind_dec,
			 ${$listref}[$i]->{'ra'},
			 ${$listref}[$i]->{'dec'});
	if ($ad<=10/3600) {
	    return $i;
	}
    }
    return -1;
}

sub bin_sort {
    my $list=shift;
    my $searchval=shift;

    my $low_index=0;
    my $high_index=$#{$list};

    my $search_low=$searchval-10;
    my $search_high=$searchval+10;
    my $done=0;
    my $high_count=0;
    my $old_high_index=0;
    my $old_low_index=0;
    while ($done==0) {
	my $mid_index=floor(($high_index-$low_index)/2+$low_index);
	if (${$list}[$mid_index]->{'ra'}>$search_high){
	    $high_index=$mid_index;
	} elsif (${$list}[$mid_index]->{'ra'}<$search_high){
	    $low_index=$mid_index;
	}
	my $diff=$high_index-$low_index;
	if (($high_index-$low_index)==2){
	    $high_index=$low_index+1;
	    $done=1;
	} elsif (($high_index-$low_index)==1){
	    $done=1;
	} elsif ($high_index==$low_index){
	    $done=1;
	}
	if (($high_index==$old_high_index)&&($low_index==$old_low_index)){
	    $high_count++;
	} else {
	    $high_count=0;
	}
	if ($high_count==100){
	    $done=1;
	}
	$old_high_index=$high_index;
	$old_low_index=$low_index;
    }
    my $real_high=$high_index;

    $low_index=0;
    $high_index=$#{$list};
    $done=0;
    my $low_count=0;
    $old_high_index=0;
    $old_low_index=0;
    while ($done==0) {
	my $mid_index=floor(($high_index-$low_index)/2+$low_index);
	if (${$list}[$mid_index]->{'ra'}>$search_low){
	    $high_index=$mid_index;
	} elsif (${$list}[$mid_index]->{'ra'}<$search_low){
	    $low_index=$mid_index;
	}
	my $diff=$high_index-$low_index;
	if (($high_index-$low_index)==2){
	    $low_index=$high_index-1;
	    $done=1;
	} elsif (($high_index-$low_index)==1){
	    $done=1;
	} elsif ($high_index==$low_index){
	    $done=1;
	}
	if (($low_index==$old_low_index)&&($high_index==$old_high_index)){
	    $low_count++;
	} else {
	    $low_count=0;
	}
	if ($low_count==100){
	    $done=1;
	}
	$old_low_index=$low_index;
	$old_high_index=$high_index;
    }
    my $real_low=$low_index;

    return ($real_low,$real_high);
}
